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The statistical properties of large, disastrous earthquakes have significance in the study of de- 
viations from the regular behavior of earthquakes. In this study, to help understand the irregular 
£SJ ■ behavior of earthquake occurrences, the large deviation function for the frequency of earthquakes is 

calculated for the first time. We study the temporal sequence of the largest earthquakes in simple 
models of earthquakes, one-dimensional forest-fire models in which the fluctuation in loading and 
C\j . fracture process are taken into consideration. We introduce four different models of fixed trigger 

sites to represent the starting points of rupture propagation. The size- frequency distributions and 
1 scaled large deviation functions for the frequency of the largest earthquakes in the system are cal- 

X^/y ■ culated and the phase tables are presented. The calculated large deviation functions are compared 

with those of the homogeneous Poisson process and of the one-site forest-fire model. We find that 
the large deviation function largely depends on the model parameters and the fixed trigger sites, and 
in most cases the large deviation function deviates from that of the homogeneous Poisson process, 
' '. suggesting that the largest earthquakes in the model do not occur randomly. The relation between 

the size- frequency distribution and the large deviation function for the frequency is discussed. 



& ; I. INTRODUCTION 

C/3 



Statistical indices for characterizing earthquakes are important for understanding the mechanism of earthquakes. 
The Gutenberg-Richter (GR) law [1] for the magnitude of an earthquake is well established in seismology, and it 
is well known that the b-value depends on both time and space. However, in the tail-part of the distribution that 
■ corresponds to large magnitudes, the functional form is still a subject of concern [2] because we have insufficient data 
from large earthquakes. Moreover, being able to estimate the frequency or probability of rare events such as disastrous 
O ■ earthquakes is a big concern from the point of view of hazard assessment. 

To further study the frequency of large earthquakes, we adopted a large deviation function (LDF) [3, 4]. The LDF 
is related to the probability of rare events, which constitute the tail part of the probability distribution. The LDF 
is universal in the sense that it is an asymptotic form whenever the number of elements is large. Recently, much 
interest has been paid to the LDF for current in nonequilibrium statistical physics, where the LDF is expected to 
OS. serve as a "thermodynamic" function [5, 6]. The LDF is also used in the analysis of activity in glassy systems, at 
the transition between active state and inactive state [7]. Besides the LDFs for current and activity, the LDF for 
frequency is also of recent concern in counting processes such as photon-counting processes [8-10]. Budini [9] studied 
the "thermodynamic" framework of a point process by using the LDF for temporal frequency. The earthquake 
sequence can be regarded as a point process such as the well-known epidemic-type aftershock sequence (ETAS) model 
[11]. The "thermodynamic" approach can be applied to such earthquake point processes. Studies that analytically 
calculated the LDF in "thermodynamics" and statistical physics, e.g., [5], preceded numerical calculations, even 
• • ■ though analytical calculation of the LDFs in complex systems is difficult. Recently, Monte Carlo methods have been 
m ^ \ introduced to obtain the LDF of a system described by master equations in discrete time [12] and in continuous time 
[13, 14]. These methods enabled us to calculate the LDFs in the complex systems of earthquakes. 

Because the recurrence time of large earthquakes is longer than several tens of years, we do not have enough data on 
large earthquakes, and observational studies of the frequency of large earthquakes are limited. In order to overcome 
this, we can generate enough data by simulating earthquakes. In this paper, we study the LDF for the frequency 
of simulated earthquakes in a one-dimensional forest-fire model, which can be understood as a minimalist model 
for earthquakes. Originally, the forest-fire model was introduced to simulate forest-fires [15]. Drossel and Schwabl 
[16] represented forest-fires in the model with four processes: planting of trees, ignition, propagation of fire, and 
extinguishing of fires. To separate the timescale of the planting process and those of the latter three, a 'lightning' 
model was introduced [17, 18]. The lightning model reduces the last three processes to a single process, the vanishing 
of a cluster of trees. 
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Besides recent applications to real forest-fires [19], forest-fire models have been applied to simulate earthquakes 
[20]. When used as an earthquake model, loading on a fault corresponds to planting a tree and triggering of an 
earthquake corresponds to igniting a hre. The idea of representing an earthquake with a randomly expanded cluster 
was introduced by Otsuka [21]. Newman and Turcotte [22] studied the cycles of large earthquakes, which were defined 
as the percolated clusters in a 2D forest-fire model. To establish a minimalist model of earthquakes and to estimate 
the predictability of earthquakes, Vazqucz-Prada et al. [23] arrived at a model similar to a ID lightning forest-fire 
model. Recently, heterogeneous configurations of ignition (trigger) sites were introduced by Tejedor et al. [24] in 
2D forest-fire models to represent the variation of faults. They found that the size-frequency distributions of the 
simulated earthquakes could be expressed by power laws. In their models, there are some asperity regions where large 
earthquakes tend to occur. According to the behavior of the size-frequency distribution, they obtained phase maps, 
which were found to be dependent on the configuration of the trigger sites. 

In this study, we adopted four ID forest-fire models with different numbers of trigger sites as simple models of 
earthquakes. We obtained the size-frequency distribution and numerically calculated the LDFs for the frequency of 
the system-size earthquakes. A system-size earthquake is an earthquake whose size is characterized by the system 
size of a model, and it is the largest earthquake in the system. We classified the behaviors of the size-frequency 
distribution and the LDFs for frequency of the system-size earthquakes into 'phases'. The LDFs for frequency of the 
system-size earthquakes were scaled by the mean frequency of the system-size earthquake. 

First, we introduce the four models used in this study and write their master equations. We give a brief introduction 
to the LDF used in this study and to the method of calculating the LDF. Next, we show the results of the size- frequency 
distributions of the four models and present a phase table. Similarly, the LDFs for the frequency of the system-size 
earthquakes scaled by the mean frequency are calculated to obtain a phase table. We discuss the relationship between 
the phase of the size-frequency distribution and that of the LDF, and we present our conclusions. 



We studied the forest-fire models in ID on a lattice of length L. To take into account the heterogeneous nature 
of faults, we introduced the four models as shown in Fig. 1, Ml, M2, M3, and MA. Ml has trigger sites only at 
the left edge, M2 can have triggers at both the edges, M3 can have triggers at both the edges and at the site m 
(2 < m < L — 1), and MA can have triggers at all sites. MA is a common homogeneous forest-fire model [25], while 
Ml, M2, and M3 may represent heterogeneous faults. In Ml, M2, and M3, earthquakes of various sizes can nucleate 
at the trigger sites, and the other sites are broken only by large earthquakes which nucleated at the trigger sites. 

To simulate a cellular automaton model such as a forest-fire model, an update rule is required [26]. Among the 
cellular-automaton-like models, the results of a random update rule are consistent with steady state solutions [27, 28] 
of the master equation in continuous time for an asymmetric simple exclusion process [29]. For the lightning model, 
a formulation of the master equation in continuous time is available [30] , and this enables us to calculate numerically 
the LDF. Thus we adopted the random update rule in the models Ml, M2, M3, and MA. 

The simulation procedure is as follows: First a site is chosen at random. If the site is empty, the site is loaded 
at the rate p per unit time. If the site is loaded and the chosen site is a trigger site, an earthquake is triggered at 
the rate / per unit time. When an earthquake is triggered, triggering also occurs at the neighboring sites until an 
unloaded site is encountered. The probabilities of the loading and the triggering processes describe the fluctuation of 
the dynamics in the loading and rupturing processes. 

The state of the jth site (j = 1, • • ■ , L) is described by the occupation number Tj, which takes 1 when the site is 
loaded and when unloaded. Loading stress on the jth site is expressed by the transition from rj = to Tj = 1, and 
stress release (an earthquake) by the transition from rj = 1 to Tj = 0, where ' represents the configuration before the 
transition. When an earthquake is triggered at one of the loaded sites from j + 1 to j + s, and the sites j and j + s + 1 
are empty, an earthquake of size s occurs, and it is expressed by the transition from {rj, • • • , rj +s+1 } = {0, 1, • • • ,1,0} 
tO {Tj,-- - ,T j+s+1 } = {0,--- ,0}. 

For simplicity, we write the configuration of the system C = {ti, • • • , tl}- We introduce a probability P(C; t) that 
the system is in the configuration C at time t. The master equation is written as 



II. MODEL AND LARGE DEVIATION FUNCTION 



A. 



One-dimensional forest-fire models 




C'^C 



(1) 



where W(CC') is the transition rate from C to C. The summation Ylc^c re P resen ts the sum over all configurations 
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FIG. 1: Description of models Ml, M2, M3, and MA with lattices of length L. Trigger sites are represented by lightning marks. 



of C except for C. W(CC') is written as 



W{CC) = p y £["^r i . 1 r } _ 1 (l-7i)T j 5 Tj+lT}+l ---]+W f {CC% 
3 = 1 



(2) 



where Wj{CC) is the trigger term that depends on the model, and S xx > is the Kronecker delta. The parts with a 
suffix less than 1 or greater than L are replaced by unity (e.g., S T0T ^ is replaced by 1.). The exact forms of Wf(CC) 
are written as 



W f (CC) = /]Tt((1 - n) • --rja - r,)(l - rj +1 )(l - r j+1 )8 Tj+2T , +2 
i=i 



for Ml, 



W f (CC') = W 2 = fJ2^-ri)---r' j (l-T j )(l-T^ 1 )(l-T j+1 )6 Tj+2 



3+2 



for M2, 



for M3, and 



for MA. 



i=i 



m L 

W f {CC>) = ^ 2 +/^^---V 2 ^ 2 [(l-rj_ 1 )(l-r,- 1 )rj(l-r,). 

j = l k—m 

T 'k( l - T k)(l - 7fc +1 )(l - Tfe+l)](5 Tfc + 2 ^ + 2 • • • 
L — 1 m L 

W>(C<7') = W 2 + /^^^...^._ 2ir ,_ 2 [(l-rj_ 1 )(l-T j _ 1 )rj(l-r j 

m— 2 j — 1 k—m 

T-'d 1 - T fe)(! ~ ^fc+iX 1 - T k+i)]5 Tk+2 T- k+2 ■ ■ ■ 



1 L-j-2 



(3) 



(4) 



(5) 



(6) 



B. Large deviation function 

The mean frequency of earthquakes of size s per unit time x(s) is written as 



x(s) 



N(s) 



(7) 
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Here, N(s) is the number of earthquakes of size s for elapsed time t. We may write N(s) as 

N(s)=J^X(s;C i+1 ,C i ), (8) 

i=l 

where X(s; Cj+i, Cj) is 1 when an earthquake of size s occurs and is otherwise, when the configuration changes from 
Ci to Ci+i, and Nt is the total number of configuration changes. The probability that the mean value x(s) is equal 
to the frequency x is written as P(x). P(x) is asymptotically written as 

P(x) ~ e-^ {x) (9) 

for large t, where the function <j){x) is called the large deviation function (LDF). The LDF satisfies <j>(x m (s)) = 0, which 
indicates the statistical mean x(s) converges to the true mean x m (s) when t — > oo, and <j>(x) > 0. For earthquakes, 
X(s; C2, Ci), • ■ • , -X"(s; Cn t1 Cn t -i) are not independent and identically distributed (i.i.d.) random variables. For 
i.i.d. random variables we have the central limit theorem. In contrast, when we study earthquakes, the LDF becomes 
significant because of the non-i.i.d. nature. 

A large deviation function <p(x) has a corresponding 'generating function' fi(X), where A is the conjugate variable 
of x. The generating function fi(X) is defined as 

e /x(A)t = ( e A*t) _ J e (>-x-4>{x))t dx . (10) 

4>{x) and fjb{\) are related by the Legendre transform as 

4>(x) = max[a;A — yu(A)]. (11) 

For example, the large deviation function for the frequency of events in a point process that obeys a homogeneous 
Poisson process is 

X 

<pp(x) = xlog(-) - x + a, (12) 
a 

and the corresponding generating function is 

fip(X) = a(e x - 1), (13) 

where a is the rate of event occurrence and the suffix P represents the Poisson process. 

To estimate the generating function /it(A), we calculate the modified master equation written as 

^-P x (C;t) = [W x (CC')P x {C';t) - W(C'C)P x (C;t)}, (14) 

where P X (C; t) satisfies the differential equation (14) under the initial condition P A (C*o; to) = P{Co',to) with an initial 
configuration Co, and W X (CC) — W{CC')e xx ^ C ' C K This modified master equation is also expressed by a matrix, 
and the largest eigenvalue of the matrix is asymptotically equal to /i(A) when t is large. We show here the calculation 
of the one-site forest-fire model to illustrate the derivation of /u(A). The transition matrix of the one-site forest-fire 
model H is written as 

H _ (-W({1}{0}) W({0}{1}) \_(-P f \ (1 v 

{ w({i}{0}) -w({o}{i}) )-{ p -f)> (lb > 



and the modified transition matrix H A as 



(16) 



l_iA _ ( -^({1}{0}) W({0}{l})e xx ^^ \_(-P fe 
~ \ W({1}{0}) -^({0}{1}) J { P -f 

By calculating the largest eigenvalue of this modified matrix, the generating function of the one-site system yUi(A) is 
given by 



Mi(A) = ^ 



-P~f+ yj(p-f) 2 +*pfe x 



(17) 
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FIG. 2: (a) Comparison of P(x) obtained by the normal Monte Carlo (M.C.) method, the cloning Monte Carlo method, and 
the matrix method, with different time steps, (b) Comparison of (j>(x) obtained by the same methods as in (a). 



The LDF is calculated from the generating function using the Legendre transform relation (11). By estimating A* that 



satisfies jfr(x\ — /x(A))|a=a» = 0, the LDF for the frequency of earthquakes of the one-site system <j>\{x) is obtained 



as 



<t>i{x) = xlog 



^-(2x +v /4.x2 + (p-/)2) 



p+f 



f) 2 + 8^ 2 + (p - /) 2 . (18) 

The mean frequency of earthquakes in the one-site system is xi m (l) — ^r* ■ This is the matrix method for obtaining 
the LDF. 

In the present study, the largest eigenvalue of the modified transition matrix is calculated numerically for L < 14 to 
obtain the generating functions. A cloning Monte Carlo method can also be used to calculate the generating function 
[13] for larger system sizes, where the matrix method cannot be used because the requirements of memory size and 
computation time would be unrealistic. However, the generating function is calculated more accurately by the matrix 
method than by the cloning Monte Carlo method. Details of the cloning Monte Carlo method used in this study are 
given in several references [13, 14, 31], so we will just give an outline of the method. For a given initial state, we 
prepare Nc copies (clones), each with its own clock. In each transition step, the earliest clone is chosen, and it evolves 
to a new configuration at a rate proportional to the modified transition rate. The amount of time that has elapsed 
from the previous transition and the number of the chosen clone would be copied or pruned (y) are calculated after 
the transition. During the copying and pruning process in this procedure, the number of clones Nc is kept constant. 
The generating function is calculated iteratively from the values Nc and y. 

The cloning Monte Carlo method and the matrix method can calculate the LDF in a wider range of x than can the 
normal Monte Carlo method, which is the method we used to simulate the forest-fire model. Values of P(x) and <fc(x) 
calculated by the three different methods are shown in Fig. 2 for the model M2 with p = 1.0, / = 0.1, and L = 12. For 
the normal Monte Carlo, 2 14 and 2 17 time steps and 2 14 ensemble members were simulated, and P(x) was obtained 
by making a histogram of x. 2 11 clones were used in the cloning Monte Carlo. P(x) and <f>(x) of the normal Monte 
Carlo with 2 14 time steps exist for 0.11 < x < 0.125, which is shown in Fig. 2, and the simulation result for the other 
values of x does not exist within the ensemble of the present calculations. This finite limiting range becomes narrower 
with increasing time steps for the normal Monte Carlo. The plotted range for the 2 17 time steps is narrower than that 
for the 2 14 time steps. Fig. 2 clearly shows that the results obtained from the three methods are in close agreement. 
This result indicates that the methods are giving sufficiently accurate values that are mutually consistent. 



III. SIMULATION RESULTS 
A. Size-frequency distributions 



In this study, the loading rate p is fixed 1.0, which corresponds to the rescaling of time by p. Fig. 3 shows the 
size-frequency distributions of simulated earthquakes for the models Ml, M2, M3, and MA, with L = 128. The trigger 
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FIG. 3: Frequency-size distribution of earthquakes for the models Ml, M2, M3, and MA with L = 128, p = 1.0, / = 1.0, 0.1, 0.01, 
and 0.001. The solid line denotes 1/s. The system-size earthquakes are highlighted by the orange dotted circles. Na is the 
total number of earthquakes. 



sites for M3 were located at sites 1,32, and 128. The trigger rate / varies as 1.0,0.1,0.01, and 0.001. We took 2 27 
Monte Carlo steps and recorded the number of earthquakes by size. The solid line in each graph of Fig. 3 denotes the 
function 1/s, where s is the earthquake size. The frequencies of the system-size earthquakes are significantly large, 
with peaks at s = L, except for MA with / = 0.1 and 1.0. 

Previous studies have found three types of behaviors ('phases') in the size-frequency distribution in the spring-block 
models and cellular automaton models of earthquakes [2, 24, 32, 33]. Following them, we classify the size-frequency 
distributions of the simulated earthquakes into three phases as follows: When the size-frequency distribution, excepting 
the system-size earthquakes, can be expressed by a power- law (the GR law), the behavior of the distribution is called 
'critical'. When the frequencies of large earthquakes whose sizes are close to L are higher than would be expected by 
the power law fit to the smaller earthquakes, it is 'supercritical'. When the frequencies of large earthquakes whose sizes 
are close to L are lower than those expected by the power law, it is 'subcritical'. These behaviors in the size- frequency 
relations of earthquakes are observed in different areas and in different faults [2], though the phase description of the 
size-frequency distribution is not commonly used by seismologists. 

The phases of the size-frequency distributions thus classified arc summarized in Table I. For Ml (Fig. 3(a)), the 
size- frequency distributions show the critical phase, independent of /. The peak of the frequency at s = L is more 
significant for smaller /. This is because all the sites tend to be loaded before triggering for small /. The exponent 
of the decay of frequency with increasing s is approximately —1.0 for / = 0.1 and 0.01. The decay rate is larger for 
larger /, because smaller earthquakes occur more easily for larger /. 

For M2 (Fig. 3(b)), the supercritical phase is observed for / = 1.0 and 0.1, where the frequency of earthquakes 
increases with s for s > 100, while the critical phase is observed for / = 0.01 and 0.001. Significant peaks at s = L 
are observed for all the values of /. For M3 (Fig. 3(c)), we observe peaks at s — 32,96, and 128 for / = 1.0,0.1, 
and 0.01. These peaks correspond to the distances between the trigger sites. The peaks at s = 32 and 96 tend to 
be unclear for smaller /, because more sites tend to be loaded before triggering and so the system-size earthquakes 
become prominent. The peaks of frequencies at s = 32 and 96 in M3, and the supercritical behavior for M2 and M3, 
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TABLE I: Phases of the size-frequency distributions in Fig. 3, where super, sub, and critical denote the supercritical, subcritical, 
and critical phases, respectively. 
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FIG. 4: Mean frequency of system-size earthquakes x m (L) versus the system size L. (a) Models Ml, M2, and M3, with 
/ = 1.0,0.1, and 0.01. (b) Model MA with / = 1.0,0.1, and 0.01. 



may be explained as follows: Small earthquakes occur frequently and release stress near the trigger sites, and clusters 
of loaded sites tend to be generated between the trigger sites. These clusters correspond to high frequencies of large 
earthquakes. 

For MA (Fig. 3(d)), the subcritical phase is observed for all cases. The size-frequency distributions may be approx- 
imated by power laws for small s. The frequency decreases more rapidly with an increase in s than would expected 
from the power law, and this rapid decrease starts at smaller values of s for larger values of /. The peak of the 
frequency of earthquakes of s — L was found for / = 0.01 and 0.001, although it was not found for / = 1.0 and 0.1. 



B. Scaled LDF of system-size earthquakes 



Next, we calculated the LDFs for frequencies of the system-size earthquakes for the models Ml, M2, M3, and 
MA, varying the system size L and the triggering rate /. Because the LDF takes the minimum value at the true 
mean frequency of the system-size earthquakes x m (L), we first evaluated x m (L) for each case of simulation. x m (L) 
depends on /, L, and the choice of model (Ml, M2, M3, or MA). Fig. 4 shows x m (L) versus the system size for 
/ = 1.0,0.1, and 0.01, and for models Ml, M2, M3, and MA. For M3, the trigger sites are located at both the ends 
and at m = L/2 + 1 for even L and m = (L+ l)/2 for odd L. Note that the results for Ml, M2, and M3 are plotted 
on logarithmic coordinates (Fig. 4(a)), while those for MA are on semilogarithmic coordinates (Fig. 4(b)). The time 
step was taken as 2 24 in each case. For Ml, M2, and M3, with / = 1.0 and 0.1, x m (L) decreased with L. When 
/ = 0.01, x m (L) seems to be independent of L, as shown in Fig. 4(a), although it approximately obeys power-law 
decay with an exponent of about —0.01 for Ml. For MA and / = 1.0, x m (L) exponentially decreases with increasing 
L. For / = 0.1 and 0.01, x m {L) takes the maximum values at L ~ 5 and L ~ 27, respectively, and exponentially 
decreases with L for large L. The decrease of x m (L) with an increase in L occurs because the preparation time for a 
system-size earthquake increases with L. 

Because x rn (L) depends on L, appropriate scaling is required to compare the LDFs for different system sizes. To 
introduce scaling by x m (L), which is simply written as x m hereinafter, we define a scaled variable z as 



(19) 
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In this scaling, the LDF of the Poisson process is given by 

^M = (z + l)log(z + l)-z. (20) 

•Em 

The division by x m on the left-hand side of (20) is unnecessary if the time t is scaled by x m , and the LDF written by 
z omits x m and can simply be written as <pp{z). The introduction of z enables us to compare the calculated LDF to 
the LDF of the Poisson process (20) more clearly. We used two methods to calculate x m , the Monte Carlo method, 
which was used for the results shown in Fig. 4, and fitting with the least-squares method. We prefer the least-squares 
method because we want the result to be self-consistent. To obtain x m by the least-squares method, for each case, 
we numerically calculated (j>(x) and fit it with a polynomial, then used the x that minimizes this polynomial as x m 
to scale the LDF. The order of the polynomial was taken to be 20. 

The scaled LDFs for the numerically calculated frequencies of system-size earthquakes are shown in Fig. 5. The 
system size L varies as L = 6, • • • ,12, 16, and 32. The LDFs in the cases of L = 16 and 32 are obtained by the cloning 
Monte Carlo method, while the others are by the matrix method. It is difficult to calculate the LDF for MA with 
/ = 1.0 and L = 16 because the number of system-size earthquakes is very small, and so instead, the LDF for L = 13 
is plotted. The red solid curve in each panel denotes the LDF of the Poisson process 4*p{z) (20), the blue dashed curve 
denotes the scaled LDF of the one-site system <f>i(z), and the black dashed curve denotes the quadratic function as 
a reference. In calculating 4>i(x), we introduced the effective triggering rate /' of the system-size earthquake, which 
is assumed to be proportional to the number of trigger sites n t , and so /' = n t f. When the LDF is a quadratic 
function, P(x) is Gaussian. It is noted that P(x) becomes Gaussian when the earthquake sequence is i.i.d. Scaling by 
x m worked well for some cases. The scaled LDFs for different values of L collapse onto a curve for Ml with / = 0.01 
(Fig. 5(c)) and for MA with / = 1.0 (Fig. 5(j)). In contrast, in the other cases, the LDFs for different values of L arc 
scattered. 

From the comparison of the obtained LDFs of the system-size earthquakes with those of the Poisson process, we 
find three phases of behavior: Poisson, half-Poisson, and non-Poisson. The LDFs for MA with / = 1.0 are well 
approximated by the LDF of the Poisson process, independent of L (Fig. 5(j)); this is called the 'Poisson phase'. In 
this case, the occupation of all sites rarely occurs because of frequent triggering. For Ml, M2, and M3, with / = 0.01, 
4>p(z) fairly approximates the simulated LDFs for z < 0, while it deviates from the LDFs for z > (Figs. 5(c), (f), 
and (i)). We call this behavior the 'half-Poisson phase'. In the half-Poisson phase, the plots are better approximated 
by <f>i(z) than by (f)p(z). (f>i(z) is close to 4>p{z) for z < 0, though it deviates from (f)p(z) for z > 0. (f>i(z) gives a 
good approximation in the half-Poisson phase. This is because earthquakes of s < L rarely occur for small /, and 
the system-size earthquakes are dominant. When the frequency of system-size earthquakes is much larger than that 
of other earthquakes, the system is approximated by the transition between two states {0, • • • ,0} and {1, ■ • • , 1}, 
which is similar to the one-site forest-fire model of two states {0} and {1}. For z < 0, the number of system-size 
earthquakes during a time interval t is small compared to the mean, and the time intervals between successive system- 
size earthquakes are longer than the average. When the time interval is long, all the sites tend to be loaded before 
triggering. This results in nearly random occurrence of system-size earthquakes, and the LDF may be approximated 
by the LDF of the Poisson process. It is remarked that the origin of the Poisson behavior in the half-Poisson phase 
for z < is different from that of the Poisson phase. The Poisson behavior in the half-Poisson phase originates in the 
Poisson process of triggering, while in the Poisson phase it originates in the random occurrence of full loading of the 
system. 

For the other cases, the simulated LDFs clearly deviate from that of the Poisson process, and this is called the 
'non-Poisson phase'. Although the LDFs show non-Poisson behavior for MA with / = 0.1 and 0.01 (Figs. 5(k) and 
(1)), we anticipate that the distributions converge to the Poisson curve with increasing L, as discussed in the next 
section. For this reason, we write N(P) for MA with / = 0.1 and 0.01. The phases of the LDFs thus classified are 
shown in Table II. For Ml, M2, and M3, we do not have a clear explanation for the non-Poisson behavior. The 
troughs in the LDFs for Ml with / = 1.0 and 0.1 (Figs. 5(a) and (b)), and for M2 and M3 with / = 0.1 (Figs. 5(e) 
and (h)), are deeper for larger L. In contrast, the troughs in the LDFs for M2 and M3 with / = 1.0 (Figs. 5(d) and 
(g)) are shallower for larger L. A deep trough in the LDF indicates that the frequency of system-size earthquakes 
is commonly close to x m , and as the trough becomes deeper, the system-size earthquakes occur more periodically 
because the earthquake sequence is a renewal process. A renewal process is characterized by a distribution of time 
intervals between successive events (here the events are the system-size earthquakes) . The sequence of the system-size 
earthquakes in the present forest-fire models is a renewal process because the system always becomes empty after each 
system-size earthquake. The LDF for the frequency of a renewal process depends on two factors: the time-interval 
distribution and the trajectory of the time intervals (see the appendix). Let us assume an LDF which has a very deep 
trough at x m , and a distribution of time intervals between successive system-size earthquakes with the dominant peak 
at At. In this case, P(x)(~ e - *"^)) is close to a delta function with the peak at x m , which means that the probability 
of the trajectories of system-size earthquakes deviating from x = x m is very small. The trajectories that give x = x m 
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FIG. 5: The plot of <j>(z) for models Ml (/ = 1.0(a), 0.1(b), 0.01(c)), M2 (/ = 1.0(d), 0.1(e), 0.01(f)), M3 (/ = 1.0(g), 0.1(h), 
0.01(i)) and MA (/ = 1.0(j), 0.1(k), 0.01(1)). 



consist of almost equal time intervals ~ At or of time intervals that contain deviations from At. The latter case is 
unlikely to occur because a deviation from At must be supplemented by other earthquakes that also deviate from At. 



IV. DISCUSSION AND CONCLUSIONS 

A. The size-frequency distribution and the LDF of system-size earthquakes 

The LDFs for MA with / = 0.1 (Fig. 5(k)) and 0.01 (Fig. 5(1)), which belong to N(P), display a clear failure of the 
scaling by the mean frequency x m . For MA with / = 1.0, (i) the LDF shows the Poisson phase (Fig. 5(j)); (ii) there 
is no peak at the system-size earthquake in the size- frequency distribution (Fig. 3 (d)); and (iii) the mean frequency 
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MA 
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N(P) 
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TABLE II: Classification of the LDF by / of the four models. P, N, and H denote the Poisson, non-Poisson, and the half-Poisson. 
The N(P) phase becomes Poisson when L is large enough. 




FIG. 6: Size-frequency distributions of MA for different L with (a) / = 1.0, (b) / = 0.1, and (c) / = 0.01. 



of system-size earthquakes is approximated by an exponential function of L (Fig. 4(b)). The characteristics (ii) and 
(iii) hold for / = 0.1 and 0.01 when the system size is large enough, as we show in detail below. 

Fig. 6 shows the size-frequency distributions of earthquakes for MA with / = 1.0,0.1, and 0.01. The distributions 
for MA with / = 1.0 do not show peaks at s = L, while those for / = 0.1 do show peaks at s = L. For MA with 
/ = 0.1 and 0.01, the peak of frequency at s = L is clear when L is small, but it becomes unclear with an increase in 
L. This suggests that the characteristic (ii) is satisfied for large L. 

Fig. 4(b) shows that the mean frequency of system-size earthquakes behaves asymptotically as an exponential 
function of L for MA with large L, indicating that the characteristic (iii) holds for large L. Thus we confirm that 
the characteristics (ii) and (iii) are satisfied for large L and MA with / = 0.1 and 0.01, and we expect that the 
characteristic (i) is also satisfied in these cases. 

In the model MA with large L and /, earthquakes are frequently triggered before the system is fully loaded. This 
results in the occurrence of many small earthquakes and few system-size earthquakes, leading to characteristics (ii) 
and (iii). The configuration of fully loaded system appears rarely and randomly, and therefore the sequence of the 
system-size earthquake follows approximately the Poisson process. Thus we expect that the LDFs of system-size 
earthquakes in the model MA may be categorized into the Poisson phase for large L. Note that this expectation has 
not been proven theoretically; it is left for a future study. 
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So far we have discussed relations among the LDF for the frequency of system-size earthquakes, the size-frequency 
distribution of earthquakes, and the dependence on system size of the mean frequency of system-size earthquakes 
for the model MA, which suggest that the characteristic time sequence of system-size earthquakes is related to the 
size-frequency distribution. What about the other models? For the non-Poisson behaviors of Ml, M2, and M3, we 
do not find clear relations between the size-frequency distributions and the LDFs for the frequency of system-size 
earthquakes. For M3 with / = 1.0, with increasing L, the LDF approaches the LDF of the Poisson process (Fig. 5(g)), 
which may be similar to the N(P) phase. However, in this case, neither the characteristic (ii) nor (iii) is satisfied, 
suggesting that the behavior is different from that of the model MA. 

B. Conclusion 

The size-frequency distributions and the LDFs for the frequency of system-size earthquakes for four kinds of ID 
forest-fire models were calculated numerically. The LDFs for the frequency of the system-size earthquakes mostly 
deviate from the LDF of the Poisson process. We classified the behaviors of the size-frequency distributions into 
three phases: supercritical, critical, and subcritical. We also classified the behaviors of the LDFs for the frequency 
of the system-size earthquakes into three phases: Poisson, non-Poisson, and half-Poisson. The Poisson phase of the 
LDF is related to the subcritical phase, where the peak at the system-size earthquake is unclear in the size-frequency 
distribution. This relation has yet to be confirmed by the calculation of the LDF for large L where the mean frequency 
of system-size earthquakes exponentially decreases with an increase in system size. The calculation of the LDF for 
such a large L is difficult at present. 

For real seismic activities, the statistical properties of the frequency of large earthquakes of magnitudes greater 
than 7.0 was studied in [34], and it was found that the sequence of large earthquakes obeys the Poisson process for 
magnitudes over 7.3. Although the ID forest-fire models are too simple to characterize the complex properties of real 
earthquakes, the present study may provide some insight into the relationship between the size-frequency distribution 
of earthquakes and the statistical properties of interevent times of system-size earthquakes. Whether the system-size 
earthquakes occur at random, obeying the Poisson process, or with some regularity, is related to fault properties. 

The present study focused on earthquakes. Deviation from the Poisson process appears in many contexts of science, 
such as bunching in the photon- counting process and a spike train of neurons. The LDF approach can also be applied 
to those phenomena, and we hope this study furthers understanding of the underlying physics of systems described 
by point processes. 
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Appendix A: LDF for the frequency of a renewal process 

The sequence of system-size earthquakes in the ID forest-fire models can be regarded as a renewal process because 
each system-size earthquake releases all the stress from the system, leading to an empty configuration. Note that the 
statistics of the occupation time of a renewal process is discussed in [35] . 

We consider the probability p(N, t) that N events occur during a time interval t. p(N, t) can be regarded as a 
function with frequency x — N/t, and the function can be written as P(x). The LDF for the frequency (f>(x) behaves 
asymptotically as the probability P(x), which is expressed by 

Mx) = lim --\ogP(x). (Al) 

t— >oo t 

We consider the sequence of interval times between successive events t\, ■ ■ ■ , tjv-i, and the probability density function 
w(U). For a renewal process, the probability of earthquakes occurring with intervals t\, ■ ■ ■ , ijv-i is written as 



N-l 

Probfa, • • • ,t;v-i) = Y\ w ( t i)- 

i=l 



(A2) 
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Since Prob(t\, ■ ■ ■ , tjv-i) is the probability for the trajectory of a sequence of earthquakes, we need to sum up the 
probabilities of all the trajectories to obtain the probability of N earthquakes in the time interval t. Thus p(N, t) is 
given by the integral with t\, ■ ■ ■ , tjv of Prob(t\, • • • , fjv) times the probability of no earthquakes during t — X^i* 
for i > X)i=i *»• is written as 

t ft N-1 N-1 N-1 

p(N, t) = / dt N _! ■■■ dti TT [«;(*<)] ^(t - V ti)0(t - V ti)), (A3) 
J o i= i j= i i= i 

where 0(t) is a step function satisfying 9(t) = 1 for t > and 9(i) = otherwise. 9(i — Y^i=i t-i) selects the 
trajectories between t > Y^a=i ti an d * < Si^i 1 VK*) is the probability of no earthquakes during t. Here we note 
that 

/oo 
di'w(*') (A4) 

is satisfied. 

From the equations (A3) and (A4), we see that p(N,t) is determined by the two factors: w{t) and the trajectories 
of the time intervals t\, ■ ■ ■ , tjv-i that satisfy t > J^iLi 1 
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